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Black hole production in tachyonic preheating 
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We present fully non-linear simulations of a self-interacting scalar field in the early universe 
undergoing tachyonic preheating. We find that density perturbations on sub-horizon scales which 
are amplified by tachyonic instability maintain long range correlations even during the succeeding 
parametric resonance, in contrast to the standard models of preheating dominated by parametric 
resonance. As a result the final spectrum exhibits memory and is not universal in shape. We find 
that throughout the subsequent era of parametric resonance the equation of state of the universe 
is almost dust-like, hence the Jeans wavelength is much smaller than the horizon scale. If our 2D 
simulations are accurate reflections of the situation in 3D, then there are wide regions of parameter 
space ruled out by over-production of black holes. It is likely however that realistic parameter values, 
consistent with COBE/WMAP normalisation, are safetly outside this black hole over-production 
region. 

PACS numbers: 98.80 Bp, 98.62 Lv, 04.70 Bw 



I. INTRODUCTION 

Over the past decade, the theory of reheating after in- 
flation has been developed and it was recognized that 
reheating can be accompanied by an era of preheating, 
where fluctuations of scalar fields are greatly amplified 
either by parametric resonance 0, 0, El 0, H| or by tachy- 
onic instability 0, S 01 ■ in preheating, energy transfer 
from the inflaton field to another field occurs in bursts 
where non-perturbative effects are important for the evo- 
lution of scalar fields. 

A common feature among the various preheating mod- 
els is that preheating leaves an imprint on the subsequent 
evolution of the universe, by amplifying the density per- 
turbation especially on sub-horizon scales at that time by 
many orders of magnitude. If the resultant density per- 
turbations at horizon crossing becomes of order unity, 
black holes can be formed by gravitational collapse |j- 
Since the energy density of produced BHs increases pro- 
portional to the scale factor relative to the total energy 
density of the universe in the subsequent radiation dom- 
inated era, BHs will dominate density of the universe at 
late time even if a tiny fraction of energy is converted 
into BHs at preheating. There are strong constraints on 
the abundance of BHs for a wide range of mass scales 
|lfj| and they can be a powerful means for constraining 
inflationary models 0| ■ 

In the past, it was pointed out by several authors that 
BHs can be over-produced during preheating for two- 
field models 0, 0, |2^|. In these models, parametric 
resonance causes the rapid amplification of field fluctua- 
tions. Then, non-perturbative effects such as rescattering 
may be crucial for calculating the abundance of produced 
BHs. Indeed, it was shown via lattice simulations that 
BHs are unlikely to be over-produced during preheating 
[l4j . Accurate simulations are required since a small er- 
ror in the estimate of backreaction time which ends the 



resonance causes a large error in the amplitude of den- 
sity perturbation at the end of preheating because field 
fluctuations grow exponentially during preheating. 

In this paper, we use lattice simulations to examine the 
evolution of density perturbations during tachyonic pre- 
heating and the succeeding parametric resonance, mainly 
focusing on the production of black holes. To evalu- 
ate the amplitude of density perturbations taking the 
non-linear interactions into account, we performed two 
and three-dimensional lattice simulations using the C++ 
code LATTICEEASY developed by Felder and Tkachev 
0. 

The equations of motion of the minimally coupled 
scalar field that we solve are 
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where p, a and H are, respectively, the total energy den- 
sity, the scale factor and the Hubble parameter, and (■ • •) 
denotes spatial average. To estimate the variance of den- 
sity perturbations at the horizon scale (k = aH), we 
simply use the power spectrum Vs(aH). Here Vg(k) is 
defined by 
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Vs(k)S(k-k'), 
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and 5k is the Fourier component of the density contrast 
6 = 6p/(p). 



II. TACHYONIC PREHEATING 

It was noted in the lattice simulations of 0, Q that 
the tachyonic instability in the context of spontaneous 
symmetry breaking can be the most efficient process of 
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energy transfer from the potential energy of scalar field 
into the energy of its inhomogeneous modes. If V"((f>) 
is negative, fluctuations with comoving momenta k < 
\/\V"\ grow exponentially and can often become quite 
large. As a result of backreaction, the spatial average of 
the field experiences only a few oscillations around the 
minimum of the potential or may not even oscillate at 
all. 

One important characteristic in tachyonic instability is 
that the longer wavelength modes have the largest growth 
rate as long as the metric perturbations can be neglected, 
which will be a good approximation for the dynamics in- 
side the horizon. On the other hand, it is well known that 
curvature perturbations on a constant energy density hy- 
persurface remain constant on super-horizon scales if the 
evolution of the energy density is described solely by a 
single equation of state. Thus no tachyonic instability oc- 
curs far outside the horizon in single scalar field models 
|20j , though we expect that perturbations on sub-horizon 
scales are amplified by tachyonic instability, which may 
still lead to over-production of BHs. Therefore it is in- 
teresting to see how large the density perturbations can 
be amplified become through the tachyonic instability on 
sub-horizon scales where one can neglect metric pertur- 
bations. 

Here we consider Coleman- Weinberg potential j2l| 

VW = \tf(tog$-\)+±to\ (4) 

which was used in the original model of new inflation 
[22I l2"il ] . Recently it was found in [25| that preheating in 
this model is so efficient that the perturbative description 
of the reheating process pjjl |29| does not apply. 

To see for which parameters of the potential BHs are 
most likely to be produced, let us consider the growth of 
the density contrast 5 in the linear approximation. Ap- 
proximately when the background 4> reaches the value 
4> = v(v/M p i), the negative mass squared of the infla- 
ton field, — V", exceeds the squared Hubble parameter, 
H 2 w Au 4 /M 2 ( . After this epoch, tachyonic instability 
sets in. Therefore, when <f> — <p , one can still set the 
initial amplitude 8<p to ~ H ~ \f\v 2 /M p i at around the 
horizon scale (k/a ~ H). On large scales the growing 
mode of S(j> grows in proportion to <f>. Neglecting the 
expansion of the universe, which will be valid when the 
mass scale exceeds the Hubble scale (<f> > <fi ), one can 
say that <f> grows in proportion to 4> 2 from energy conser- 
vation j25j. Therefore at a small k we have 8<f> cx 4> 2 . 

Using the relations mentioned above, one can deduce 
8(f) /<f) ~ \^XM pt (j)/v 2 . In the parameter region \f\ < 
(v/Mpi), S(j)/(j) stays less than unity. Therefore fluctu- 
ations do not become non-linear. In this case the am- 
plitude of density perturbations also will not reach the 
0(1) values required for BH production. On the other 
hand, in the case where s/X > (v/M p i), fluctuations of 
the inflaton field evolve into non-linear regime. 

If we assume that primordial density perturbations ob- 
served, e.g. in the anisotropy of the cosmic microwave 



background, are produced from the vacuum fluctuations 
of 0, A is fixed to about 1CT 12 from COBE (COsmic 
Background Explorer) normalization. If we adopt this 
value, the most interesting parameter for v is 10~ e M p i 
or less from the above argument. But in this case, 
there is a large hierarchy (six orders of magnitude) be- 
tween the minimal wave number for tachyonic instabil- 
ity to occur, which is determined by the horizon scale 
(~ y/Xv(v/M p i)), and the one for parametric resonance 
to occur (~ \f\v). Therefore we cannot incorporate these 
two scales simultaneously into our numerical simulations 
for the parameter range, v J$ yXM p i, in which large am- 
plification of fluctuations is expected via tachyonic insta- 
bility. Hence instead of doing simulations for such a small 
value of v/Mpi, we adopt A = 10~ 4 and v = I CP 2 M p i, 
so that both two scales are included in our simulations. 
We note that in the double inflation scenario j23, 
A > 10~ 12 can be realized. 

We start our simulations at the time when d> = v-S— . 

T M p i 

The initial condition for Vs^ is determined by solving the 
linearized equations, 

6<f> k + 3H^ k +(^ + V"(<f>)\ S<j> k = 0. (5) 

from <f> <C 0o t° 4> — 4>o- We mainly present the re- 
sults for two-dimensional simulations (2048 2 ) with trans- 
lation symmetry along the z-axis imposed to cover the 
still widely separated typical length scales of the differ- 
ent instabilities discussed above. Later on, we briefly 
comment on three dimensional cases. 

III. NUMERICAL RESULTS 

A. A = 10- 4 , v/M p i = 10~ 2 

This is a marginal case in which non-linear interactions 
become efficient just when the inflaton reached the min- 
imum of the potential. But as we will see, simulations 
show that the linear regime lasts a little bit longer until 
the middle stage of the parametric resonance. Fig.n an d 
Fig. |21 respectively, show the evolution of Vs<f> /v and Vg 
for A = ICP 4 and v/M p i = 10~ 2 . Here the proper time 
t is measured in units of the inverse of the typical mass 
scale of the inflaton potential = s/\v. Let us first 
focus on Fig. ^ We find that long wavelength perturba- 
tions are selectively amplified during the tachyonic phase 
as expected from the linear analysis. After 4> starts to 
oscillate around the potential minimum (r > 90), para- 
metric resonance occurs at scales around k/a ~ m^. A 
remarkable thing is that once parametric resonance com- 
mences, the amplitude of Stp with long wavelengths de- 
creases in time, and its spectrum eventually obeys the 
power law cx k 2 . This scaling can be interpreted as a 
result of losing correlations on large scales by parametric 
resonance on small scales. 

Next we focus on fig. [21 We find that Vg around hori- 
zon scales (k ~ aH) becomes as large as 10 -3 , which 
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FIG. 1: Time evolution of Vs^/v for the Coleman- Weinberg 
potential for A = 10~ 4 and v/M p i = 10~ 2 . r is a dimension- 
less time defined by r = t/m^, where t is cosmological time. 
The plot (also Fig. 2) is made by averaging ten different sim- 
ulations in order to reduce the statistical errors. 
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FIG. 2: Time evolution of Vs for the Coleman- Weinberg po- 
tential for A = 10 -4 and v/M p i = 10 -2 . The horizon scale 
is located in the range —1.8 < log 10 fc/m^, < —2.0 (slightly 
outside the range of our simulation) during preheating. 



is in good agreement with the prediction of the linear 
approximation. In the linear approximation, we have 
5 ~ -pH^N ~ SN on large scales, where SN is the pertur- 
bation of the e-folding number on flat hypersurfaces and 
is roughly estimated as y/\ at the initial time {<f> = <fi a ). 
Here we used the fact that SN is conserved on large scales 
for single field models. Although the linear perturbation 
theory breaks down after the tachyonic instability, lattice 
simulations suggest that there is no further amplification 
of S on sub-horizon scales by non-linear effects in this 
case. 

Immediately after the parametric resonance starts, the 
universe begins to expand almost like the dust-dominated 
universe because the energy density of the universe is 
dominated by the long wavelength modes compared with 
the mass scale, which are selectively excited by the para- 
metric resonance. Figure shows the evolution of scale 
factor a and w respectively, where w = P/ p. As shown 
in [l7| by the extrapolation of the linear approximation, 
the lowest smoothed density contrast 5 C above which BHs 
are formed is about w. Numerical simulations [Tsl Il9l | 



which studied BH formation in the radiation dominated 
era showed 0.3 < 5 C < 0.5, supporting the conclusions of 
[l7j . In the following estimate of the produced BH abun- 
dance, we adopt this threshold so that we assume the 
relation S c = w. Assuming that the probability distri- 
bution of density perturbations at the Jeans wavelength 
~ y/wH~ l is Gaussian, the observational upper limit of 
the standard deviation of this distribution a is about 
0.1S C for a wide range of BH mass scales. We find from 
Fig. 13 that the time average of w is smaller than 10~ 2 . 
Therefore the threshold value on Vs is ~ (0.1u>) 2 10~ 6 , 
which is much smaller than the values we find from our 
numerical simulations; see fig. In this case, BHs are 
over-produced. 




FIG. 3: Evolution of the scale factor a(r) and w. 
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FIG. 4: Comparison of two-dimensional simulations with dif- 
ferent grid sizes (256 2 ) and (2048 2 ). In both cases the time 
evolution of Vs is plotted for A = 10 -4 and v/M p i = 10 -2 . 

We also performed simulations for different values of 
v/M p i around 10~ 2 such as 5x 10~ 2 , 5x 10~ 3 , 2x 10~ 3 and 
so on. All the cases exhibit similar behavior as Figs. 
and □ 

During the tachyonic phase, 5 also starts to grow fol- 
lowing the amplification of 5(f) at large scales, and the 
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spectrum Vs becomes slightly red. Interestingly, at the 
late phase of the parametric resonance (r > 120 in 
fig. EJ fluctuations of the 0-ficld on large scales decrease 
and tend to obey the scaling Vs<p oc fc 2 , while density 
perturbations on large scales remain close to the flat 
spectrum. After r > 145, Vs becomes almost time- 
independent. This will be explained by causality. Energy 
transfer faster than the velocity of light cannot occur. If 
we neglect the presence of correlations between different 
Fourier modes of 5(f>, the red spectrum seen in Vs on large 
scales at a late time after parametric resonance cannot be 
explained from the blue spectrum of Vs^jv This means 
that there are strong correlations among different Fourier 
modes. As is expected, if we increase the amplitude of 
initial fluctuations on large scales by hand, the resulting 
Vs also becomes larger. 

We also performed three-dimensional simulations with 
grid size (256 3 ), which is not sufficient to include both 
the scales of tachyonic instability and of parametric res- 
onance. Before showing the result of three-dimensional 
simulations, we present comparison of two-dimensional 
simulations with different grid sizes (256 2 ) and (2048 2 ) 
in figQ] In the former case the scales at which paramet- 
ric resonance is efficient, k/a ~ m^, are not fully covered 
by the simulation. Then, the amplitude on large scales is 
relatively enhanced (but only by a small factor), presum- 
ably because the tachyonic preheating is kept efficient for 
a longer period due to weaker backreaction. Fig.[S]shows 
Vs after the tachyonic instability for three-dimensional 
case. We see that this spectrum is almost identical to 
the two-dimensional case with TV = 256 presented in 
hg- El and its growth has already been saturated on large 
scales. This suggests that three dimensional simulations 
with 2048 3 should yield a spectrum similar to the corre- 
sponding two-dimensional one and thus there is a strong 
potential of over-producing BHs in the three-dimensional 
case. 




k/m, 

FIG. 5: Three-dimensional time evolution of Vs for A = 10 
and v/Mpi = 10~ 2 . 



B. A = 1, v/Mpi — 1CF 2 

Fig.[|)]and[3show the evolution of Vs^/v and Vs respec- 
tively for A = l,v/M p i = 10~ 2 . In this case, non-linear 
interactions are already important during the tachyonic 
phase. We have removed the contribution to (cj> 2 ) from 
the vacuum fluctuations of the short wavelength modes 
by introducing the cutoff, Vs<p — ► Vs<f>/(1 + (k/k c ) a ). 
This cutoff is necessary because otherwise the dynam- 
ics is largely affected by the shortest wavelength modes 
which are not resolved well numerically 

From fig. we see that, contrary to the case with 
A = 10~ 4 and v/M p i = 10 -2 where linear approximation 
is valid during the tachyonic phase and the subsequent 
parametric resonance occurs efficiently, parametric reso- 
nance is less efficient in this case because non-linear in- 
teractions among modes around k < have become 
important when inflaton reaches the potential minimum. 
In particular, Vs^/^ik) at small k does not approach the 
scaling regime given by Vscf, oc fc 2 , and long range corre- 
lations are still maintained after parametric resonance. 

Next we focus on the evolution of the density pertur- 
bations. We found that Vs at a small fc almost remains 
constant after r ^ 90. The value of Vs at small fc is about 
0.1 which is about one order of magnitude smaller than 
the one in the linear approximation. This shows that the 
early importance of non-linear interactions suppress the 
terminal amplitude of density perturbations presumably 
due to the effect of backreaction rather than increase it. 

Before closing this section, we briefly comment on the 
possibility of producing BHs for other parameter regions 
of A and v/M p i which are not covered in the simulations. 
The results of the numerical simulations of both A and 
B show that the non-linear interactions which may be- 
come important during the phase of the tachyonic insta- 
bility or that of the parametric resonance do not amplify 
the density perturbations and hence the final curvature 
perturbations on constant energy density hypersurface 
do not exceed the one evaluated by the linear perturba- 
tion theory where the curvature perturbations is given 
by - V\. 

If we assume that this result holds for smaller values 
of A and v/M p i which are not covered in the simulations, 
then the amplitudes of perturbations evaluated by the 
linear perturbation theory, which is v'A in this case, will 
determine the abundance of produced BHs. If we also 
assume that the threshold value of over-production of 
BHs docs not depend on A and v/M p i, then BHs will 
indeed be over-produced for A > 10 -6 . 



IV. DISCUSSIONS 

In summary, we studied non-linear evolution of den- 
sity perturbations on sub-horizon scales by two and 
three dimensional fully nonlinear lattice simulations. We 
adopted the Coleman- Weinberg potential as a model of 
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FIG. 6: Time evolution of P S(f> / v for A = 1 and v/M pi = 10" 2 . 




FIG. 7: Vs for A = 1 and v/M p i = 10~ 2 . 



tachyonic instability, since this model is expected to cap- 
ture typical features in more extended preheating models. 

We found that the density perturbations below the 
horizon scale are selectively amplified by tachyonic in- 
stability and thereafter maintain their amplitude, while 
fluctuations of the inflaton field decrease in time if the 
following parametric resonance dominates the perturba- 
tions. If the non-linear interactions become significant 
during the tachyonic phase, parametric resonance is less 
efficient and long range correlations of inflaton fluctua- 
tions do not decay. 

We also found that there is no enhancement of per- 
turbations due to the non-linear interactions and the 
terminal value of the spectrum Vs on subhorizon scales 



is bounded from above by the linear estimate SN 2 ~ 
A. Hence if we extrapolate this result to the realis- 
tic A ~ 1CP 10 determined by temperature fluctuations 
of the cosmic microwave background, we see that BH 
over-production is unlikely to occur for realistic A with 
SN ~ 1CT 5 . 

We also found that the resulting power spectrum of 
density perturbations strongly depends on the initial con- 
ditions when preheating is dominated by the tachyonic 
instability. The latter point is in contrast to the cases 
dominated by parametric resonance, in which the spec- 
trum tends to obey a universal power law Vs oc k D (D is 
the spatial dimension) independently of the initial shape 
of the spectrum. 

Finally we briefly comment on the role of metric per- 
turbations which are not taken into account in our sim- 
ulations. It is necessary to include metric perturbations 
for the evolution of long wavelength perturbations. In 
linear perturbation theory, the evolution equations for 
inflaton fluctuations on flat slicing is given by |30j 

8<t> + 3HS6 + K.S6 + V"S6 = a' 3 ( —<j>A 56. (6) 
a 1 \H J 

The right hand side of eq. iJSJ represents the contributions 
from metric perturbations. We solved eq. JBJ) numerically 
and compared the resulting power spectrum with the one 
without metric perturbations. We found little difference 
between the two power spectra, which shows that the ne- 
glecting the metric perturbation is a good approximation 
within linear perturbation theory. 

There is a possibility that non-linearity may enhance 
the perturbations when metric perturbations are taken 
into account but this possibility is beyond the scope of 
the current work. 
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